// tspcluster - Traveling Salesman Problem
// Copyright (C) 2006  Frederik Carlier
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

using System;
using System.Text;
using System.Diagnostics;
using System.Collections.Generic;
using System.Collections.ObjectModel;

namespace UGent.Tsp
{
    [Serializable()]
    public class Results
    {
        #region Constants
        // By default, disable the histogram
        public const int DefaultHistogramNumberOfBins = 0;
        public const int DefaultCorrelationCount = 500;
        #endregion

        #region Constructors
        public Results()
            : this(DefaultCorrelationCount)
        {
        }

        public Results(int correlationCount)
            : this(correlationCount, DefaultHistogramNumberOfBins)
        {
        }

        public Results(int correlationCount, int histogramNumberOfBins)
            : this(correlationCount, histogramNumberOfBins, 0.0, (double)histogramNumberOfBins)
        {
        }

        public Results(int correlationCount, int histogramNumberOfBins, double histogramLower, double histogramUpper)
            : this(correlationCount, histogramNumberOfBins, histogramLower, false, histogramUpper, false)
        {
        }

        public Results(int correlationCount, int histogramNumberOfBins, double histogramLower, bool histogramLowerExtended, double histogramUpper, bool histogramUpperExtended)
        {
            // Initialize correlation
            this.correlationCount = correlationCount;
            if (correlationCount > 0)
            {
                this.correlationProduct = new double[correlationCount];
                this.buffer = new double[correlationCount];
                this.firstValues = new double[correlationCount];
            }

            this.transitions = new Transitions(histogramNumberOfBins, histogramLower, histogramLowerExtended, histogramUpper, histogramUpperExtended);
        }
        #endregion

        #region Fields
        // If two results were merged, some things don't work
        bool merged;

        // General statistics
        double sumOfValues;
        double sumOfSquaredValues;
        double minimum = double.MaxValue;
        double maximum = double.MinValue;
        double current; // The last added value
        int currentBinIndex = -1;

        // The number of measurements
        // that are in the statistical sample;
        long count;

        // The number of rejected measurements
        long rejected;

        // Correlation
        [NonSerialized()]
        double[] firstValues;
        [NonSerialized()]
        double[] correlationProduct;
        [NonSerialized()]
        double[] buffer;
        [NonSerialized()]
        int correlationCount;

        // Histogram and TransitionMatrix
        Transitions transitions;


        #endregion

        #region Properties
        // The percentage of moves that were accepted
        public double AcceptanceRate
        {
            get
            {
                return (double)(count - rejected) / (double)count;
            }
        }

        public int CorrelationCount
        {
            get
            {
                if (merged)
                    throw new InvalidOperationException();

                return correlationCount;
            }
        }

        public long Count
        {
            get { return count; }
        }

        public double Minimum
        {
            get
            {
                return minimum;
            }
        }

        public double Maximum
        {
            get
            {
                return maximum;
            }
        }

        public Transitions Transitions
        {
            get
            {
                return transitions;
            }
        }
        #endregion

        #region Methods
        public void Merge(Results results)
        {
            if (results == null)
                throw new ArgumentNullException("results");

            this.sumOfValues += results.sumOfValues;
            this.sumOfSquaredValues += results.sumOfSquaredValues;
            this.minimum = Math.Min(this.minimum, results.minimum);
            this.maximum = Math.Max(this.maximum, results.maximum);
            this.count += results.count;
            this.rejected += results.rejected;
            this.transitions.Normalize();
            this.transitions.Merge(results.transitions);

            // Other values don't work.
            merged = true;
        }

#warning This API needs refactoring
//        [Obsolete("Use Add(double, TransitionType) instead.")]
        public void Add(double value)
        {
            Add(value, TransitionType.Success);
        }

        // A rejected value still has valuable information about the transition probabilities
        // starting at the current energy.
//        [Obsolete("Use Add(double, TransitionType) instead.")]
        public void Add(double current, double rejected)
        {
            Add(rejected, TransitionType.Failed);
        }

        public void Add(double value, TransitionType type)
        {
            if (count == 0 && type == TransitionType.Failed)
                throw new InvalidOperationException();

            // Set the "new value", to add to the statistics. If the transition is
            // succesful, this should be "value". If failed, get the previous value --
            // stored in "current"
            double newValue = type == TransitionType.Success ? value : current;

            sumOfValues += newValue;
            sumOfSquaredValues += newValue * newValue;
            current = newValue;

            count++;
            if (type == TransitionType.Failed)
                rejected++;

            // The index of the current value
            long index = count - 1;

            if (newValue > maximum)
                maximum = newValue;

            if (newValue < minimum)
                minimum = newValue;

            if (correlationCount > 0)
            {
                // Save the value in the buffer for the correlation values;
                long bufferIndex = index % correlationCount;
                buffer[bufferIndex] = newValue;

                // Save in the first values buffer if required
                if (index < correlationCount)
                    firstValues[index] = newValue;

                // Calculate the product of the i'th and i+k'th values, required
                // for calculating the autocorrelation.
                long dUpper = Math.Min(count, correlationCount);

                // Normally, we do this
                // inside the loop:
                // int k = Cyclic.GetIndex(index - d, correlationCount);
                // but it uses quite some CPU cycles. Hence,
                // we optimize. K starts as index. After every loop;
                // k goes to k-1. If k < zero, we have to increment.
                int k = Cyclic.GetIndex(index - dUpper, correlationCount);
                for (int d = 0; d < dUpper; d++)
                {
                    // Get the index of count - d = index - d
                    correlationProduct[d] += buffer[bufferIndex] * buffer[k];
                    k--;
                    if (k < 0)
                        k += correlationCount;
                }
            }

            if (transitions.Dimension > 0)
            {
                int binIndex = transitions.GetBinIndex(value);

                if (count > 1)
                {
                    transitions.AddTransition(currentBinIndex, binIndex, type);
                }

                if (type == TransitionType.Success)
                {
                    currentBinIndex = binIndex;
                }
            }
        }

        public double GetAverage()
        {
            return sumOfValues / (double)count;
        }

        // Gets the deviation of the sample population. Is an estimator for the
        // total population deviation, sigma.
        public double GetDeviation()
        {
            if (count < 2)
                return 0;

            // Dit levert de steekproefvariantie
            // http://nl.wikipedia.org/wiki/Variantie
            //   \sqrt{\frac{1}{N-1} \left( \sum_{i=1}^{N} x_i^2 - N \avg{x}^2 \right)
            // = \sqrt{\frac{1}{N-1} \left( \sum_{i=1}^{N} x_i^2 - \frac{1}{N} ( \sum_{i=1}^{N} x_i )^2 \right)
            return Math.Sqrt((sumOfSquaredValues - sumOfValues * sumOfValues / (double)count) / ((double)(count - 1)));
        }

        // Gets the mean square error on the average, assuming a standard distribution of
        // the population. 
        // http://en.wikipedia.org/wiki/Mean_squared_error
        public double GetMeanSquareError()
        {
            return GetDeviation() / (double)Count;
        }

        public double GetRootMeanSquareError()
        {
            return Math.Sqrt(GetMeanSquareError());
        }

        public double GetAutocorrelation(int d)
        {
            if (merged)
                throw new InvalidOperationException();

            double average = GetAverage();
            double deviation = GetDeviation();

            double numerator, denominator;

            double sumOfFirstValues = 0;
            double sumOfLastValues = 0;
            for (int i = 0; i < d; i++)
            {
                sumOfFirstValues += firstValues[i];
                sumOfLastValues += buffer[Cyclic.GetIndex(count - i - 1, correlationCount)];
            }

            numerator = correlationProduct[d] - average * (2 * sumOfValues - sumOfFirstValues - sumOfLastValues) + (double)(count - d) * average * average;
            denominator = (double)(count - 1) * deviation * deviation;

            return numerator / denominator;
        }

        public double GetAutocorrelationTime()
        {
            if (merged)
                throw new InvalidOperationException();

            // rho(x) = exp(-x/tau)
            // ln(rho(x)) = -x/tau
            // tau = -x/ln(rho(x))
            //
            // If the autocorrelation time is 0 (uncorrelated series),
            // the previous calculation won't work.
            // 
            // Also, for large values of x (d), the autocorrelation time
            // is approximately zero and we are left with statistical noise.
            // Hence, once we encounter a correlation with absolute value
            // smaller than 0.01, we stop the calculation.


            int max = 100;

            for (int d = 1; d <= max; d++)
            {
                if (Math.Abs(GetAutocorrelation(d)) < 0.01)
                {
                    max = d - 1;
                    continue;
                }
            }

            if (max == 0)
                return 0;

            double[] x = new double[max];
            double[] y = new double[max];

            for (int d = 0; d < max; d++)
            {
                x[d] = (double)(d + 1);
                y[d] = -1.0 * (double)(d + 1) / Math.Log(Math.Abs(GetAutocorrelation(d + 1)));
            }

            // The values y[d] should be constant, so we can do a zero-th degree polynomal fit.
            double[] coefficients = LeastSquares.Fit(x, y, 0);

            return coefficients[0];
        }
        
        public Results Clone()
        {
            Results value = new Results(this.correlationCount);
            value.buffer = this.buffer == null ? null : (double[])this.buffer.Clone();
            value.correlationProduct = (double[])this.correlationProduct;
            value.count = this.count;
            value.current = this.current;
            value.currentBinIndex = this.currentBinIndex;
            value.firstValues = this.firstValues == null ? null : (double[])this.firstValues;
            value.maximum = this.maximum;
            value.merged = this.merged;
            value.minimum = this.minimum;
            value.rejected = this.rejected;
            value.sumOfSquaredValues = this.sumOfSquaredValues;
            value.sumOfValues = this.sumOfValues;
            value.transitions = this.transitions.Clone();
            
            return value;
        }
        #endregion

        #region Static Methods
        public static double[,] ToArray(Results[,] value)
        {
            if (value == null)
                return null;

            int length0 = value.GetLength(0);
            int length1 = value.GetLength(1);

            double[,] output = new double[length0, length1];

            for (int i0 = 0; i0 < length0; i0++)
                for (int i1 = 0; i1 < length1; i1++)
                    output[i0, i1] = value[i0, i1] == null ? 0.0 : value[i0, i1].GetAverage();

            return output;
        }
        #endregion
    }
}
